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Abstract 

We study the phase diagram of the site-diluted Ising model in a wide di- 
lution range, through Monte Carlo simulations and Finite-Size Scaling tech- 
niques. Our results for the critical exponents and universal cumulants turn 
out to be dilution-independent, but only after a proper infinite volume ex- 
trapolation, taking into account the leading corrections-to-scaling terms. 

PACS: 75.50.Lk. 05.50.+q, 68.35.Rh, 75.40. Cx, 
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1 Introduction 



The magnetic phase diagram and critical properties of many magnetic materials 
can be described by a simple model (the Heisenberg Hamiltonian): 

H= E W S l (!) 

Here Sf is a spin operator (Latin indices refer to lattice sites, while Greek ones 
represent the spin components). J^g is a coupling matrix which is usually short- 
ranged and can represent either the Ising, XY or Heisenberg models by properly 
dealing with the spin index. One can understand eq. (Q) on the basis of the 
exchange interaction between the electrons of the external shells of the atoms. This 
interaction, expressed as in eq. (Q), is symmetric in the spin index. Nonetheless 
if one puts the atoms on a crystalline lattice, the material tends to magnetize in 
the so-called axes or planes of easy magnetization given by the symmetry of the 
crystal. 

One typical example are the uniaxial crystals, as the hexagonal lattices, where 
the magnetization can choose as subspace of easy magnetization the c axis or its 
orthogonal plane. In the first case the system is well described assuming that the 
magnetic momenta point in the c direction and it should be described by the Ising 
model. In the second one, the material should be studied by means of the XY 
model. 

However no pure material exists in Nature, and it is mandatory to consider the 
effects of non-magnetic impurities. The simplest way to do so, is by considering a 
modified version of ([!]) 

H = E J%ZitiS?sP, (2) 

where the e's are quenched, uncorrelated random variables, chosen to be 1 with 
probability p (the spin concentration) or with probability 1 — p (the impurity 
concentration, or spin dilution). The rationale for the quenched approximation is 
that usual relaxation times for the non-magnetic impurities are much longer than 
the corresponding for spin dynamics. For non-frustrated systems, the phase dia- 
gram of (|2|) in the temperature-dilution plane consists of a magnetically disordered 
(paramagnetic) region at high temperature, separated from an ordered (ferromag- 
netic) region at lower temperatures. The dilution dependent critical temperature, 
T c (p), obviously equals the pure model value at p = 1. It lowers for larger dilution 
values, until the extreme case T c (p c ) = at the site percolation threshold for the 
concentration of the magnetic atoms. 

Not many general results have been obtained for the Hamiltonian (Q). The 
most popular one is doubtless the Harris criterion [1]. This criterion states that 
the critical behavior of (||) will be the same as for (I) if the specific- heat critical 
exponent, a, is negative, while a new Universality class will appear if a > 0. In 
the latter case its is possible to show jz] that a for the diluted model is negative: 
one of the effects of the dilution is to smooth the critical behavior of the system. 
The only one between the generic models for magnetism (Ising, XY, Heisenberg) 
displaying a > in three dimensions, is the Ising model. 

There are other physical contexts in which the Hamiltonian (||) has been stud- 
ied. For instance, its four dimensional Ising version has been recently investigated 
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(see m and references therein) in connection with the puzzling problem of finding 
non asymptotically- free interacting theories in four dimensions. The two dimen- 
sional model is also interesting as a playground for exactly solvable field-theories, 
and has also been considered (see [Q, ||, |(| and references therein) . 

As already stated, the materials displaying Ising-like behavior in very pure 
samples should behave differently when the impurities concentration increases. In 
fact, according to Harris, an infinitesimal impurity concentration should be enough 
to spoil the Ising behavior. However this will happen in very narrow intervals of 
temperatures, which may be unreachable experimentally. 

The Hamiltonian (0) can be studied in the low dilution regime by means of 
analytical perturbative renormalization-group methods JtL |[ E|. They find a new 
fixed-point, thus implying that the critical exponents along the T c (p) line are 
dilution independent and different from their pure Ising value. The predicted 
correlation length exponent, v, ranges from 0.697 in ref. to 0.67 in ref. [[)). 
This should be compared with the Ising model result, v = 0.6300(15), given in 
ref. For the order parameter critical exponent, (3, a value of 0.35 is obtained, 
contrasting with 0.3250(15) for pure Ising ]ic|j . For the susceptibility exponent, 
the perturbative analysis predicts 7 = 1.32 and the corresponding Ising one [jl0| is 
1.241(2). This is maybe the quantity more easily comparable with experimental 
results, as the magnetic susceptibility can be very precisely measured. 

The study of the Hamiltonian (|^) beyond the low disorder regime, is restricted 
to the Monte Carlo (MC) method. Many simulations have been performed in the 
last seventeen years 0, [TJ, |l|, |ll[ [0| [l6| . The first study, on small lattices |ll| 
was compatible with the new fixed point scenario. However further simulations 
[ fi2| found results rather suggesting a continuously varying value of the critical 
exponents along the critical line. A Monte Carlo Renormalization Group study Jl4|] 
found a value for the v exponent consistent with the perturbative one at p = 0.8. 
However, for p — 0.9 their results did not differ from the pure Ising model, while 
for p < 0.8 they could not find meaningful results. More recent simulations 13]] 
suggested a single fixed point scenario with v = 0.77(4), confirmed in ref. 1 15f| 
where v — 0.78(1) was found at p = 0.4. This puzzle of mutually contradicting 
results started to make sense in ref |l^] . In this work, the crucial observation that 
the exponents measured in a finite lattice are transitory was made. Unfortunately 
the statistical errors at large dilution did not allow for a definite conclusion. 

When writing this paper a new MC work on this model has appeared |l7j . They 
obtain v = 0.682(2) at p — 0.8 but a markedly different result (y = 0.717(8)) at 
p = 0.6. 

In this paper we present the first sound numerical evidence for a random fixed 
point dominant along the whole critical line. This is achieved by means of a Finite- 
Size Scaling (FSS) analysis, in a wide dilution range (0.4 < p < 0.9, the percolation 
threshold being at p c w 0.31 |l8|]). The investigation of very diluted samples is 
made possible by a p-reweighting method, which allows to extrapolate the simu- 
lation results obtained at p to a close p' value jl9[ ||, [I] . A careful consideration 
of the scaling corrections is needed, in order to get the right value in the infinite 
volume limit. In this system, the first corrections-to-scaling exponent, u>, is very 
small (lu « 0.4, see ref. (?]]). Thus, the confusing results in previous MC studies 
can be understood as an unusually large contribution of the scaling corrections. 
After a proper consideration of this problem, we find dilution independent critical 
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exponents in quantitative agreement with perturbative calculations. 

Other theoretical problem of interest is the absence of self-averaging at the 
critical point. This means that the disorder-realization variance of quantities such 
as the magnetic susceptibility or the specific-heat, at the critical point, is a fixed, 
non-zero fraction of their mean values even in the thermodynamical limit. It 
has been argued pC| that this fixed fraction is an universal number. In ref. 
this fraction for the susceptibility is calculated analytically and numerically in four 
dimensions. In this work, we numerically calculate this ratio, along the critical line 
T c (p). After the compulsory infinite volume extrapolation, an universal, dilution 
independent result is found. A very recent simulation ||l7f has questioned the 
universality of these ratios. However, these authors do not perform any infinite 
volume extrapolation, making their conclusions necessarily not definitive. 

The experimental study is still not completed. For instance, indications of 
the expected new universality class were obtained in the Ising antifcrromagnet 



Fei_a;Zn x F2, studied in the reduced temperature range 10 3 < t < 10 1 |2l[ 22 



In this system the order parameter exponent was found to be [3 — 0.36 [21 
while the obtained susceptibility exponent, 7, was 1.44(6)|2^]. Also in ref. p2] 
cusp-like behavior of the specific-heat was found, so no divergence was expected. 
This yields v > 2/3 through standard hyperscaling relations. Another system in- 
vestigated was a dysprosium aluminum garnet doped with yttrium ]23[ ] , for which 
(3 = 0.385(25) was obtained at a 5% dilution. The results regarding the (3 exponent 
have been questioned in ref. [ p4| where Mno.5Zno.5F2 was studied by synchrotron 
magnetic X-rays scattering. These authors conclude that the experimental errors 
to date are too big to distinguish between the pure Ising and the diluted (3 values. 
Maybe the strongest evidence found for a new Universality Class has been re- 
ported in ref. [ p5[ studying Mni_ x Zn x F2 by means of neutron scattering. Critical 
exponents v = 0.70(2) and 7 — 1.37(4) were found. 

The layout of the paper is as follows. In section 2 we define the model and the 
observables to be measured in the numerical simulation. In section 3 we provide the 
necessary technical details about the MC methods. Section 4 is devoted to Finite 
Size Scaling techniques. After that, in section 5, we present our numerical results 
and discuss the need for an infinite-volume extrapolation. This is considered in 
section 6. We present our conclusions in section 7. 



2 The Model 

We have considered the site-diluted Ising model on the single-cubic lattice, with 
nearest neighbors interaction. We will work in a lattice of linear size L, with 
periodic boundary conditions. The Hamiltonian is 

H = -(3 22 eiCjCiVj , (3) 

<i,j> 

where a are the usual Z2 spin variables. The e's are the quenched random variables 
introduced in (j|). We shall refer to an actual {e^} configuration as a sample. We 
study the so-called quenched disorder: that is, for every observable it is understood 
that we first calculate the average on the {<7i} variables with the Boltzmann weight 
given by exp(— if), the results on the different samples being later averaged. 
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To avoid confusions, wc will denote the Ising average with brackets, while the 
subsequent sample average will be overlined. The observables will be denoted with 
calligraphic letters, i.e. O, and with italics the double average O = (O). The total 
nearest-neighbor energy is defined as 

£ = e ^ e i a i ■ ( 4 ) 

(id) 

The energy is extensively used for extrapolating the results obtained for an observ- 
able, O, at coupling /3 to a nearby 0' coupling |2(j] and for calculating /3-derivatives 
through its connected correlation. For instance, one can define the specific-heat 

as 

C = d p TT^ = - ((£ 2 )-{£) 2 ) , (5) 

V being the total number of sites in the lattice, L 3 . 
The normalized magnetization is 

i 

In terms of the magnetization we can give a convenient definition of the suscepti- 
bility as 

X = V(M 2 ), (7) 

its Binder parameter being 



94 = t;- ■ (8) 

Another kind of cumulant, meaningless for the pure system, can be defined as 



_ (M 2 ) 2 (M 2 ) 

92 — 2 ' V y J 

(M 2 ) 

This quantity would be zero in the thermodynamical limit if self-averaging is to 
be found. A very useful definition of the correlation length in a finite lattice, 
reads \M 



*H ,r a ,-,:J - do) 



X/F-I 
,4sin 2 (7r/L) 

where F is defined in terms of the Fourier transform of the magnetization 



<?(fc) = -^E eifcr£ ^ - (ii) 



as 

V 



F = — (|£(27r/£, 0,0)| 2 + permutations) . (12) 
This definition is very well behaved for the FSS method we employ [B8fl. 
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3 The Monte Carlo Update 



The method of choice for an Ising model simulation is a cluster- method p3| . The 
most efficient variety for the pure model is the Wolff single-cluster update [ [30| . 
However, in diluted systems, very small (even individual) groups of nearly (or 
completely) isolated spins can appear. These groups are scarcely changed with a 
single cluster method. Thus we have constructed our elementary MC step (EMCS) 
as 250 cluster flips complemented with a Metropolis step. For the largest dilutions 
(p = 0.4, 0.5) the presence of isolated intermediate-sized groups of spins makes the 
thermalization too slow. For these dilutions, in the EMCS we have carried out 
a standard Swendsen-Wang sweep every 200 single-cluster flips. We discard fOO 
EMCS for equilibration, then measuring after every EMCS. The autocorrelation 
times for all observables are very small (near f EMCS in the largest lattice), but we 
have also controlled that our update method correctly thermalizes, by comparing 
hot and cold starts. 

A disordered model simulation gets characterized by two parameters, the num- 
ber of samples generated (Ns), and the number of independent measures taken 
in each sample (Ni). Previous works (for instance 14 [f5|, |l6)) have chosen the 
Nj S> Ns regime. However (see fl), the optimal regime is 
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MS)'- (13) 

where aj is the mean variance in a sample of the observable under consideration, 
while as is the variance between different samples. Moreover, the non-vanishing 
value of <?2 shows that the susceptibility is not a self-averaging quantity, thus 
making very dangerous the small Ns regime. In this work we have fixed Nj = 200 
and N s = 20000. For p = 0.9 we performed N s = 10000. 

In addition to the usual (3 extrapolation Q, in some cases it is useful to 
perform a p extrapolation. It can be done as we know the precise distribution of 
the densities of the actual configurations (binomial distribution). Details of the 
method can be found in ref. [|| for the same model in four dimensions. 

We remark that the large number of samples used, combined with the relative 
small number of measures, makes the (3 extrapolations biased. A proper statistical 
procedure allows to cancel the bias. We address to ref. (|] for details about the 
method we follow. 



4 Finite Size Scaling Methods 

A very efficient way of measuring critical exponents p{| follows from this form of 
the FSS Ansatz 

0{L,[3,p) = L*°l v (Fo{t{L,[3,p)/L) + 0{L-»)) , (14) 

where a critical behavior t~ x ° is expected for the operator O and Fo is a (smooth) 
scaling function. From a Renormalization Group point of view, u> is the eigenvalue 
corresponding to the leading irrelevant operator. It is very important that, in the 
above equation, only quantities measurable on a finite lattice appear. Notice that 
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terms of order £^=oo are dropped from cq. (|l4|), so we assume that we are deep 
within the scaling region. 

To eliminate the unknown scaling function, we measure the quotient 

Q =0(sL,p,p)/0(L,p,p) , (15) 

at the coupling value for which the correlation length in units of the lattice size is 
the same for both lattices. So we get 

Qo\ Qi = s = s Xo/u + 0(L-«) . (16) 

Given the strong statistical correlation between Qo and Q^, the above quotient can 
be obtained with great accuracy (in fact, in our opinion, this is the best method 
available to measure the usually tiny three-dimensional 77 exponents pj|). 

In many cases (high precision computations or small lattices), it is useful to 
parameterize the leading corrections-to-scaling, thus we need to consider in the 
analysis a behavior like 

Qo\ Q(=s =s x °/» + A%L-'» + ... . (17) 

Here the dots stand for higher-order corrections, while A® is a dilution-dependent 
slope. 

The most convenient observables to measure the two independent critical ex- 
ponents, 77 and v, are found to be 

dp£, -> x = v + 1 , 
X -> x = 1/(2-77). 



5 Numerical results 

The phase diagram of the model (|^) is shown in figure [l]. In this work we have 
simulated lattices L = 8, 16, 32, 64 and 128, at dilutions p = 0.9, 0.8, 0.65, 0.5 and 
0.4. Our procedure has been the following. For p = 0.9, 0.8, 0.65, we have chosen 
a [3 coupling value where the relation 

gL,(3, P ) e(2L,A P ) 

- L~ = 2L (18) 

approximately holds. Then, we have relied on standard reweighting methods, 
which allow to extrapolate the simulation results at coupling j3 to a close /3', 
to precisely fulfill the matching condition For very diluted systems, the 

transition line is almost horizontal (see figure |l ) thus it is more convenient to use 
a p-reweighting method to extrapolate the simulation results to a nearby p' value 
(see refs. [|l9[ ||, ||)- Therefore, we have first located the (3 values for which eq. ( |l8| ) 
holds at p = 0.4,0.5, then we have fixed this (3 value, and changed p later on. In 
this way, the true critical dilutions for fixed /3, p c (ft), differ from 0.4 and 0.5 (in 
less than a 2%). Nevertheless, we shall keep referring to them as p = 0.4,0.5 in 
tables and graphics, for the sake of clarity. 

In table | we present the results for exponents v and 77 and cumulants (74 and 
(72, using eq. (|l6|) (neglecting scaling corrections). Beware that consecutive data in 
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Figure 1: Phase diagram of the model (|3|), in the inverse temperature-dilution 
plane. The dots correspond to the simulated points, while the arrow signals the 
percolation limit (/3 = oo). 



the table are anticorrelated (the results of lattice L are used once in the numerator 
and another time in the denominator in eq. ([l6])). For the error computation we 
have used a jack-knife method with 50 bins, ensuring a 10% of uncertainty in the 
error bars. Thus, we display 2 digits in these bars if the first one is smaller than 
5. 

Notice that the exponent r\ and the cumulant 34 are, before any infinite volume 
extrapolation, quite dilution independent. This can be understood because they 
show a very mild evolution with the lattice size. On the contrary, exponent v and 
cumulant (72 show a larger dependence on the lattice size and so, an infinite volume 
extrapolation is needed before one can extract definite conclusions. Nevertheless, 
one can already guess from the table that v is surely different from the pure Ising 
value and the (72 cumulant is different from zero (there is not self-averaging). The 
latter was also observed in the same model in four dimensions H , where we found 
mean field results plus logarithmic corrections. 

Another quantity of interest is the specific-heat. As stated in the introduction, 
a is negative and no divergences are expected. This is a quite difficult behavior 
to study, because FSS investigations in other models displaying a < 0, show that 
the specific-heat at the critical point is a growing, though bounded, quantity [p8| . 
For this reason we choose to study 

AC(L) = [C(2L) - C{L)] Qi=2 . 

This quantity diverges if a > 0, tends to zero if a < and goes to a constant value 
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Table 1: Critical quantities obtained from pairs (L,2L) using eq. (|l6| ) for all the 
dilutions simulated. 





L 


p = 0.9 


p = 0.8 


p = 0.65 


p = 0.5 


p = 0.4 


V 


8 


.0171(7) 


.0219(7) 


.0284(10) 


.0296(24) 


.0322(29) 




16 


.0277(7) 


.0308(7) 


.0330(8) 


.0345(19) 


.0297(16) 




32 


.0320(9) 


.0335(8) 


.0329(9) 


.0313(11) 


.0315(17) 




64 


.0349(9) 


.0346(8) 


.0335(8) 


.0329(14) 


.0326(13) 


V 


8 


.6663(14) 


.6877(11) 


.7172(16) 


.7447(24) 


.7718(32) 




16 


.6643(14) 


.6849(12) 


.7107(18) 


.7328(22) 


.7534(32) 




32 


.6631(15) 


.6836(12) 


.7048(20) 


.7189(24) 


.7382(27) 




64 


.6644(15) 


.6864(14) 


.6996(20) 


.7118(21) 


.7182(26) 


92 


8 


.0832(10) 


.1546(16) 


.2310(25) 


.2784(24) 


.3043(24) 




16 


.0861(12) 


.1500(14) 


.2077(15) 


.2371(20) 


.2551(22) 




32 


.0918(13) 


.1474(17) 


.1920(20) 


.2138(22) 


.2296(25) 




64 


.0974(17) 


.1477(12) 


.1842(19) 


.1994(21) 


.2106(16) 


94 


8 


.7049(14) 


.6900(17) 


.6814(23) 


.6900(20) 


.6989(20) 




16 


.6926(17) 


.6818(15) 


.6809(16) 


.6871(18) 


.6958(21) 




32 


.6876(19) 


.6819(16) 


.6832(20) 


.6879(17) 


.6889(20) 




64 


.6821(16) 


.6771(18) 


.6780(17) 


.6825(19) 


.6857(22) 



if the specific-heat diverges logarithmically (a = 0). In addition, the (usually 
large) background term of the specific-heat disappears. It will be convenient to 
recall that deriving the FSS Ansatz from the Renormalization Group pl| , one 
finds a behavior for the specific-heat as L 2yT ~ d (where yx = 1/f). Therefore one 
should expect the fulfillment of hyperscaling relations for the transient exponents, 
a(L) and v{V). In figure || we plot the AC(L) values obtained. As a contrast we 
also plot the corresponding values for the pure Ising model which grow, as they 
should (the data are taken from ref. p2|). We find a decreasing value of AC(L) 
for p < 0.8, as expected. Notice that the (transient) v « 2/3 found for p — 0.9 in 
table [l], implies a = through hyperscaling relations. This is very nicely shown in 
the plot, where a constant value of AC(L,p = 0.9) is seen. Plotting AC(L) against 
L a / V would be useless, because the scaling corrections go approximately as L~ 0A , 
that is, their lattice size evolution is much faster than that of the asymptotic term. 

6 Infinite volume extrapolation 

As shown in the previous section, with our statistical accuracy the values for the 
critical exponents are seen to depend on the lattice size, so an infinite volume 
extrapolation is required (see eq. ([1^)). However, one has to decide when the dots 
in eq. ( p~7|) can be neglected. Our criterium will be the following. We perform the 
fit for lattice sizes not smaller than a given L m i n . If the fit quality is reasonable 
(i.e. a not too large x 2 /d.o.f. calculated with the full covariance matrix), we 
repeat it for lattices not smaller than 2L m i n . If this last fit is also reasonable and 
the extrapolated values are compatible in both fits, we keep the central value from 
the L m in fit, but quote error bars from the 2L mm one. 
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Therefore we need an estimate for uj. We shall obtain it from the lattice size 
evolution of the scaling functions: 





= (0 


Q 5 =2 










= .94°° + 


Q 5 =2 


= .92°° + 



9±U=2 = g? + A?L-" + ..., (19) 



A^L~ U + .... 

Then we shall use this ui value to extrapolate the critical exponents v and r\. 
A reasonable value of % 2 /d.o.f. in these fits will be a consistency condition. A 
technical point of interest is that the single universality-class scenario requires the 
infinite volume extrapolation for s x °/ 1 ' to be dilution- independent. Therefore, we 
can include data of different dilutions and lattice sizes in the fit. 

In fig. |3[ we plot the minimum of % 2 /d.o.f. in a fit to cq. ([l9|), as a function 
of lu. Several points become clear. It is obvious that is not useful at all in 
order to fix w (this is not surprising as it shows almost no scaling corrections, 
A^ 4 0). We see that including the p = 0.9 data yields an untenable fit with 
L m i n = 16. Moreover, when we study the extrapolation for Qd f3 t; = 2 1+1// ", we find 
an awful result. This could have been anticipated from figure ^, where a clearly 
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Figure 3: Minimum of \ 2 as a function of u>, for the fits of eq. (19). We also plot 
with a dashed line the corresponding quantity for the Qdp£ fit. 



non-asymptotic value for the specific-heat at p = 0.9 is seen. On the contrary, 
discarding the p — 0.9 data, reasonable fits are obtained. Thus, we conclude that 
the p — 0.9 system is still crossing-over from the pure Ising fixed point to the 
diluted one, even for lattices as large as L — 128. Finally, it is evident from the 
plot that the determination of lu can be greatly improved by means of a joint fit 
of the (72 and £/ L scaling functions. The results for this fit are shown in table ||. 
According to our, conservative, dots-neglecting criterium, we find 



lu = 0.37(6). 



(20) 



Notice that the value obtained in |Q, lu — 0.42 (without error estimation), using 
the scaling-field method for momentum-space RG equations, is compatible with 
ours. 

In table [j] we present the infinite volume extrapolation for v and 77 critical 
exponents and the #4 cumulant. We see that L m i n = 16 fulfills our dots-neglecting 



criterium for 54 and r\. 
are 



For v, L„ 



8 is found to be enough. Our final values 



v = 0.6837(24)(29), 
r) = 0.0374(36)(9), 
34 = 0.673(7)(2), 



(21) 



where the first error is statistical while the second is due to the uncertainty in lu. 
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Table 2: Results of the infinite volume extrapolation of 172 and including 
data from L > L m ; n , at p = 0.4,0.5,0.65 and 0.8. Q(x 2 , d.o.f.) is the probability 
of getting a larger x 2 m the fit. 



^min 


X Vd.o.f. 


Q 


id 


Z/L 


92 


8 


46.2/21 


0.0012 


0.430(15) 


0.5890(17) 


0.1458(17) 


1G 


15.0/13 


0.31 


0.37(2) 


0.598(4) 


0.145(3) 


32 


1.95/5 


0.86 


0.38(6) 


0.587(7) 


0.150(7) 



Table 3: Infinite volume extrapolation and fit qualities for the critical exponents, 
including data from L > L m i n , at p = 0.4,0.5,0.65 and 0.8 using eq. (|l7|). The 
second error is due to the indetermination in lu — 0.37(6). 



-^min 


Extrapolation 


XVd.o.f. 


Q 


V 8 


0.6837(10)(29) 


14.0/11 


0.24 


16 


0.6838(24)(33) 


6.26/7 


0.51 


32 


0.687(6)(2) 


4.14/3 


0.25 


V 8 


0.0419(8)(20) 


96.4/11 


< io- Lb 


16 


0.0374(12)(9) 


8.92/7 


0.26 


32 


0.0374(36)(8) 


0.18/3 


0.98 


34 8 


0.6726(21)(25) 


31.5/11 


.0001 


16 


0.6734(28)(21) 


7.95/7 


0.34 


32 


0.665(7)(3) 


1.08/3 


0.78 



From ( ^l|) we obtain 

a = -0.051(7)(9), 

P = 0.3546(18)(10), (22) 
7 = 1.342(5)(5). 

For the computation of the statistical error in [3 and 7 we take into account that 
the statistical correlation between v and 77 has turned out to be negligible. 

In figure |] we show Qq^ as a function of lo for all the dilutions. We also plot 
the corresponding values for the pure Ising model. The solid lines correspond to 
the joint fit for L m j n = 8 using the data from p < 0.8. Notice that the data are 
strongly anticorrelated, therefore the apparent \ 2 on the plot is larger that the real 
one, computed with the full covariance matrix. An analogous fit for 52 is shown 
in figure |[ We remark that the p = 0.9 data point to a maybe too low value. 
This is another signature of the crossover from the Ising fixed point ((72 = 0) to 
the diluted one. 

It is interesting to compare the values for and gi with those obtained in four 
dimensions ||. 

.94 = 0.32455, 
g 2 = 0.31024. 
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Figure 4: Qd /3 (; — 2 1+1 ^ u for the different dilutions. The solid lines correspond 
to a fit enforced to yield the same infinite volume extrapolation for p < 0.8. The 
smallest lattice in the fit is L = 8 and we use u> = 0.37. The Ising data have been 
taken from fl32]. 



Finally, we can compute the infinite volume critical couplings by studying the 
crossing points of scaling functions (as ^/L and (74) measured in lattices of sizes L 
and sL. Let A/3,f, Ap^ be the deviation of these crossing point from the infinite- 
volume critical couplings. The expected scaling behavior is |3j| : 

A#,A J £«i=^L-'*- 1 /»\ (23) 

In table || we present the crossing points of £/L and 34 for the (L,2L) pair for 
all the dilutions simulated. We find again that an infinite volume extrapolation is 
needed in order to extract the critical couplings. 



Using eq. (23) for s = 2 we perform a joint fit for both scaling functions 
(74 and £/L. For this fit we take cu + \jv = 1.83(6). The final results for the 
different dilutions studied are shown in table where two values for L mul are 
used. Let us remark that our critical couplings are compatible with the results 
in 0] (f3P=°- s = 0.28578(4), /3P= - 9 = 0.24933(3)). But we definitely do not agree 
with the value (3 P=0 ^ = 0.2857609(4) quoted in ref. p7| . This is not surprising as 
in this work the corrections-to-scaling are not considered. 



13 



0.0 

o.o 



0.1 



0.2 



L 



0.3 

— u 



OA 



0.5 



Figure 5: Cumulant 172 as a function of L~ u . The solid lines correspond to a fit 
enforced to yield the same infinite volume extrapolation for p < 0.8. The smallest 
lattice in the fit is L = 16 and we use u> = 0.37. 

7 Conclusions 

We have shown, beyond the low-disorder limit, that the diluted Ising model is in 
the basin of attraction of a single fixed point. Therefore, if randomness is to be 
modelized with eq. @ , the critical exponents of an Ising system are not those of the 
pure Ising model, but those of the random fixed-point (although this may be fairly 
hard to show in a very pure sample). To establish this result we have simulated 
in a very wide dilution range, finding a consistent picture only after an infinite 
volume extrapolation. The p = 0.9 data seem, however, to be still crossing-over 
from the pure Ising fixed-point to the diluted one in lattices as large as L = 128. 

We obtain the values of the critical exponents and universal cumulants elimi- 
nating the systematic errors coming from the leading corrections-to-scaling terms. 
The previous computations did not consider these terms and were not able to con- 
trol the corresponding systematic effects. Incidentally, most of the computations 
have been carried out at p = 0.8 as in this case the scaling corrections are very 
small, and the results in small lattices seem stable. However, even in this case the 
lack of an extrapolation produces an underestimation of the errors. 

The (dilution-independent) critical exponents are shown to be in good agree- 
ment with the series estimates R ^J, jjj. The corrections-to-scaling exponent, u>, 
is measured with a 16% error and is found to be in quantitative agreement with 
the perturbative estimate Q. The smallness of this exponent explains why this 
problem is so hard to attack numerically. In fact, the total computer time devoted 
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Table 4: Crossing points of scaling functions £/L and 34 for pairs L and 2L for 
the different dilutions. 



L 


p = 0.9 


p = 0.8 


p = 0.65 


p = 0.5 


p = 0.4 


9a 8 
16 
32 
64 


.249583(30) 
.249340(15) 
.2492901(13) 
.2492924(15) 


.286002(48) 
.285765(18) 
.285758(7) 
.2857417(25) 


.37025(13) 
.370185(36) 
.370208(16) 
.3701649(48) 


.49996(27) 
.49949(6) 
.499485(30) 
.499409(11) 


.39577(28) 
.39512(8) 
.394895(33) 
.394840(13) 


HL 8 
16 
32 
64 


.249299(26) 
.249291(12) 
.2492957(44) 
.2492901(13) 


.285690(49) 
.285708(15) 
.285745(6) 
.2857394(23) 


.36961(10) 
.369986(31) 
.370147(13) 
.3701540(44) 


.49814(17) 
.49896(5) 
.499326(21) 
.499374(9) 


.39302(19) 
.39441(6) 
.394694(23) 
.394785(10) 



Table 5: Infinite volume critical couplings estimations for the studied dilutions. 
The first error bar corresponds to the statistical fit error, the second one (almost 
negligible) is due to the uncertainty in lu + \ jv exponent. For this table we use 
w + l/i/= 1.83(6). 



-^min 


rVd.o.f. 


Pc 


/3c 


16 


0.11/3 


0.394816(H)(2) 


0.852 


32 


0.04/1 


0.394821(22)(7) 


0.852 


16 


2.93/3 


0.499413(9) (1) 


0.543 


32 


0.78/1 


0.499394(17)(4) 


0.543 


16 


5.27/3 


0.65 


0.370166(5)(1) 


32 


1.53/1 


0.65 


0.370156(8)(0) 


16 


5.41/3 


0.8 


0.2857421(30)(0) 


32 


0.27/1 


0.8 


0.2857368(47)(5) 


16 


8.45/3 


0.9 


0.2492905(19)(0) 


32 


0.03/1 


0.9 


0.2492880(30)(5) 



to this work has been about 5 Intel Pentium-Pro years. As we had already shown 
in four dimensions diluted Ising models are found not to be self-averaging at 
criticality in three dimensions (see ref. fll7| for an independent verification in three 
dimensions). This is proved by showing that the quotient between the sample- 
variance of the susceptibility and its mean-value, tends in the thermodynamic 
limit to a non-zero constant independent of the dilution (it is a renormalization- 
group invariant). This quotient is measured with a 4% accuracy after the infinite 
volume extrapolation. 
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